function [sacloc, sac_thresh, peakout, thresh_data] = EL1k_sacdetect(t, speed, MULT, minthresh, peakthresh)
%
% [sacloc sac_thresh thresh_data] = EL1k_sacdetect(time, speed, [MULT: 2.5], [minthresh: -1], [peakthresh: 60])
% time = time points
% speed = speed during trial (TRY A 3 POINT FILTFILT)
% MULT = std multiplyer (DEFAULT = 2)
% minthresh = force a fixed threshold (not suggested) -1 for auto
% peakthresh = min peak speed for a saccade 60-90 range.
%
% sac_thresh = threshold (for specified range: tr)
% thresh_data(tr)=[ # saccades start_loc(tr) end_loc(tr) ];
% if the data send in needs broken up this data is useful: # of saccades|
% indexes of where the sac_thres effected data began & end
%
% coe@queensu.ca
if nargin<2
    help(mfilename)
    return
end
if nargin<4
    MULT=2;
end
if nargin<5
    minthresh=-1;
end
if nargin<6
    peakthresh=60; %peak speed must be at least 60dps
end
resolution=min(diff(t));
%EL1k data forces recording to stop between trials so let's do this section
% by section if need be.
%temp=unique(diff(t)); should return something like [2 200+ 200+ 200+ ...]; this is the inter trial dead space.
end_loc=[find(diff(t)>50); length(t)];
starloc=[1;(find(diff(t)>50)+1)];
if length(end_loc)~=length(end_loc)
    error(['length(end_loc)~=length(end_loc) in ' mfilename]);
end
fprintf('section %3d of %3d (%5.1f%%)\n',0,length(starloc),0)
sacloc=[];
peakout=[];
sac_thresh=zeros(size(starloc));
for tr= 1:length(starloc)
    fprintf([repmat('\b',1,28) 'section %3d of %3d (%5.1f%%)\n'],tr,length(starloc), tr/length(starloc)*100)
    tspeed=speed(starloc(tr):end_loc(tr));
    fix_noise=tspeed(tspeed<50 & tspeed>0); % most people use a fixed 50dps for saccade onset. here i remove all real speeds over 60dps ...
    if minthresh >= 0
        sac_thresh(tr)=minthresh;
    else
        sac_thresh(tr)=mean(fix_noise)+std(fix_noise)*MULT;
        %sac_thres=max([20 mean(fix_noise)+(MULT*std(fix_noise))]); % and find
        %the mean and 3std to create a variable threshold to compensate for different noise levels
    end
    
    s=[0;tspeed(1:end-1);0];
    d=diff(s>=sac_thresh(tr)); % speed goes above threshold
    e=[find(d==1) find(d==-1)-1];
    e=[find(d==1)-1 find(d==-1)];% ONE POINT OUT SIDE THRESHOLD
    l=e(:,2)-e(:,1);
    e(l<(10/resolution),:)=[];% saccades must last longer than 10 ms
    %     figure(555);clf;
    %     plot(speed,'.-g')
    %     hline(sac_thresh,'k')
    %     vline(e(:,1),'b')
    %     vline(e(:,2),'r')    
%     if ~all(tspeed(e(2:end,1))>=sac_thresh(tr) & tspeed(e(2:end,1)-1)<sac_thresh(tr) )
%         warning(mfilename,' some saccade start locations are flawed')
%     end
%     if ~all(tspeed(e(1:end-1,2)+1)<sac_thresh(tr) & tspeed(e(1:end-1,2))>=sac_thresh(tr) )
%         warning(mfilename,'some saccade end locations are flawed')
%     end
    n2=size(e,1);
    if n2 > 0
        peak=zeros(n2,1);
        for i=1:n2
            peak(i)=max(s(e(i,1):e(i,2)));
        end
        fsac= peak>peakthresh; % faster than peakthresh
        if ~isempty(fsac)
            sacloc=[sacloc; e(fsac,:)+starloc(tr)-1]; %#ok<AGROW>
            peakout=[peakout; peak(fsac)];%#ok<AGROW>
        end
    else
        sac_thresh(tr)=-9999;
    end %if n2 > 0
    thresh_data(tr,:)=[ size(sacloc,1) starloc(tr) end_loc(tr) ]; %#ok<AGROW>
end % for tr= 1:length(starloc)



